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AFIT/GOR/ENS/ENC/95M-06 

ABSTRACT 

This  thesis  examined  the  use  of  response  surface  methodology  (RSM)  to  estimate  the 
parameters  of  a  finite-element  groundwater  model.  An  existing  two-dimensional,  steady- 
state  flow  model  of  a  fractured  carbonate  groundwater  system  in  southwestern  Ohio 
served  as  the  calibration  target  data  set.  A  Plackett-Burman  screening  design  showed  that 
only  four  of  the  ten  hydraulic  conductivity  zones  significantly  contributed  to  the  output  of 
the  finite-element  model.  Also,  the  effective  porosity  parameter  did  not  significantly  affect 
the  model’s  output.  Using  only  the  four  significant  hydraulic  conductivity  parameters; 
four  two-level,  four-factor  designed  experiments  were  conducted  to  exploit  the  first-order 
response  surface  defined  by  a  residual  sum  of  squares  function.  Additionally,  a  central 
composite  design  and  ridge  analysis  were  used  to  adjust  the  four  parameters  and  finally 
arrive  at  a  calibrated  model  in  a  grand  total  of  146  runs.  The  final  calibrated  model,  which 
had  an  average  head  elevation  of 292  meters,  matched  the  calibration  target  data  set  with  a 
mean  absolute  error  of  only  7  mm  over  all  524  nodes  of  the  model.  RSM  provided  an 
effective  calibration  technique  to  estimate  groundwater  flow  parameters. 
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Groundwater  Model  Parameter  Estimation 
Using  Response  Surface  Methodology 


I.  Introduction 


Background 

Water  quality  is  a  topic  of  concern  in  the  United  States  and  around  the  world. 
Fresh  water  occurs  at  the  surface  of  the  earth  in  ice,  lakes,  and  rivers  and  as  groundwater 
in  the  openings  in  rocks  and  soils.  While  surface  sources  of  water  satisfy  important 
irrigation  and  population  needs,  the  predominant  source  of  fresh  water  is  provided  by 
groundwater  systems.  To  assess  the  flow  and  quality  of  groundwater  in  a  particular 
region,  groundwater  hydrologists  are  often  asked  to  predict  the  behavior  of  groundwater 
systems.  Hydrologists  use  groundwater  models  to  predict  how  the  groundwater  will  flow 
in  a  particular  area.  These  models  are  often  numerical  models  implemented  on  computers 
to  simulate  groundwater  flow. 

These  simulation  models  are  mathematical  descriptions  of  the  groundwater  flow  in 
an  underground  system.  They  incorporate  the  relevant  physical  laws  which  control  the 
system,  as  well  as  specific  hydraulic  and  geologic  (hydrogeologic)  properties  unique  to  the 
particular  groundwater  system  of  interest.  The  hydrologist  modeling  the  system  must 
properly  adjust  the  parameters  representing  the  local  hydrogeological  properties  to  ensure 
the  model  accurately  reflects  the  flow  characteristics  of  the  system  under  study. 
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The  process  of  adjusting  the  model’s  parameters  until  its  output  (typically 
hydraulic  heads)  conforms  to  the  observed  field  heads  is  termed  calibration. 
Mathematically,  this  process  of  parameter  estimation  is  known  as  solving  the  inverse 
problem.  Methods  to  solve  the  inverse  problem  are  numerous,  and  only  a  few  have  been 
implemented  to  calibrate  groundwater  models  automatically.  In  fact,  most  practitioners 
still  use  manual,  trial-and-error  calibration  methods  (Anderson  and  Woessner,  1992:235). 

One  approach  that  has  not  been  discussed  in  the  groundwater  literature  on  model 
calibration  is  response  surface  methodology  (RSM).  RSM  incorporates  a  number  of 
statistical  techniques  used  in  the  empirical  study  of  the  relationships  between  one  or  more 
output  responses  and  a  group  of  input  variables.  By  systematic,  simultaneous  adjustment 
of  multiple  hydrogeological  parameters,  the  RSM  technique  should  yield  a  calibrated 
groundwater  model  with  known  statistical  properties. 

Problem  Statement 

The  goal  of  this  research  was  to  demonstrate  that  RSM  could  serve  as  an  effective 
calibration  technique  to  estimate  groundwater  flow  parameters. 

Research  Objectives 

To  accomplish  the  stated  goal,  the  following  objectives  were  established: 

1 .  Obtain  a  data  set  of  hydraulic  heads  to  serve  as  a  calibration  target  data  set.  The  data 
set  used  by  Smith  and  Ritzi  (1993)  was  chosen  as  the  calibration  target. 
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2.  Ensure  the  calibration  target  set  produces  results  similar  to  those  obtained  by  Smith 
and  Ritzi  when  their  simulation  model  is  implemented  on  the  AFIT  computer  system. 

3 .  Conduct  a  screening  analysis  to  determine  the  statistically  significant  parameters  for 
this  groundwater  model. 

4.  Using  only  those  parameters  found  to  be  statistically  significant,  demonstrate  that 
response  surface  methods  can  be  used  to  calibrate  the  groundwater  model. 

Scope  and  Limitations 

1 .  This  study  used  the  entire  stratigraphic  section  of  the  selected  aquifer  to  demonstrate 
the  applicability  of  RSM  to  calibrating  a  groundwater  model  in  a  realistic  situation. 

2.  A  data-rich  calibration  target  data  set  of  524  nodes  was  used  to  evaluate  the  match 
between  the  groundwater  model’s  output  and  the  calibration  target  set. 

3.  Only  the  groundwater  flow  model  was  calibrated  (i.e.  the  calibration  of  heads).  The 
issue  of  solute  transport  calibration  was  deferred  to  a  later  study. 
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II.  Literature  Review  of  the  Inverse  Problem 
in  Groundwater  Systems 


Introduction 

Many  groundwater  flow  models  simulate  the  flow  of  groundwater  by  computing 
the  flow  of  water  through  relatively  small  portions  (elements)  of  the  groundwater  system. 
The  results  of  the  water  flow  through  each  of  these  elements  are  combined  to  produce  a 
simulated  description  of  the  entire  groundwater  system.  This  type  of  simulation  is  called  a 
finite-element  model.  Finite-element  models  require  estimates  of  the  hydrogeologic 
parameters  to  compute  groundwater  flow  or  to  determine  contaminant  concentrations. 

Mathematically,  the  process  of  determining  the  model  parameters  is  known  as 
solving  the  inverse  problem.  “The  inverse  problem  has  been  defined  as  that  of  estimating 
the  model  parameters  from  measurements  of  the  output.  Therefore,  the  inverse  problem 
can  be  equated  to  parameter  estimation”  (Carrera,  1987:  551).  The  difficulty  of 
estimating  these  parameters  lies  in  the  mathematical  construction  of  the  model.  A  well- 
posed  model  will  have  a  unique  and  continuous  solution.  Groundwater  models  of  the  real- 
world  do  not  contain  enough  information  to  make  a  well-posed  mathematical  model. 
Finite-element,  groundwater  models  are  thus  very  difficult  to  solve  due  to  uncertainties  in 
estimating  these  parameters.  Although  the  professional  literature  contains  numerous 
methods  on  how  to  solve  the  inverse  problem.  Response  Surface  Methodology  (RSM)  has 
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not  been  applied  to  its  solution.  The  goal  of  this  thesis  was  to  demonstrate  the  usefulness 
of  RSM  in  parameter  estimation. 

This  chapter  begins  with  a  brief  overview  of  the  mathematical  concepts  associated 
with  ill-posed  mathematical  models.  The  solution  techniques  applied  to  determining  the 
parameters  of  groundwater  flow  models  are  then  reviewed  based  upon  the  solution 
approach  —  direct  or  indirect.  Finally,  the  latest  development  in  the  inverse  solution  of 
groundwater  systems,  solving  the  coupled  (simultaneous)  flow-contaminant  problem,  is 
presented. 


The  Ill-Posed  Inverse  Problem 

For  an  inverse  problem  to  be  well  posed,  it  must  meet  the  following  three 
conditions:  1)  the  parameters  are  identifiable,  2)  the  solution  is  unique,  and  3)  the  solution 
is  stable  (Carrera  and  Neuman,  1986b:21 1).  A  parameter  is  said  to  be  identifiable  if  the 
model  output  is  sensitive  to  it.  If  a  parameter  is  nonidentifiable,  then  regardless  of  the 
value  assigned  to  the  parameter,  the  model  will  produce  the  same  output.  Even  if  all  the 
parameters  are  identifiable,  the  solution  to  the  inverse  problem  may  be  nonunique. 
Nonunique  solutions  are  usually  associated  with  multiple  local  optima.  Unstable  solutions 
have  spatially  oscillating  parameters  and  depend  heavily  on  the  initial  parameter  settings. 

If  any  one  of  the  three  conditions  is  not  met,  the  problem  is  said  to  be  ill-posed. 

Groundwater  models  are  often  ill-posed.  In  groundwater  modeling,  the  known 
information  is  limited  to  the  data  gathered  from  a  limited  sampling  of  observation  wells. 
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These  wells  are  spread  out  geographically  and  the  samples  are  spread  out  over  time.  This 
data  also  has  various  uncertainties  and  errors  associated  with  the  data  collection  process. 
These  limitations  in  the  observed  data  can  result  in  non-unique  or  unstable  solutions. 
Additionally,  when  the  hydrologist  selects  parameters,  a  thorough  sensitivity  analysis  or 
screening  should  be  performed  to  determine  which  of  these  parameters  are  identifiable  or 
significant.  However,  a  screening  for  significant  parameters  is  often  not  performed 
(Carrera  and  Neuman,  1986b:212-217).  Thus,  the  limitations  of  the  available  data, 
combined  with  the  failure  to  screen  for  the  significant  hydrologic  parameters,  render  many 
real-world  groundwater  inverse  problems  ill-posed. 


Parameter  Estimation  Methods 

In  an  attempt  to  deal  with  the  difficulties  associated  with  ill-posed  models, 
numerous  mathematical  methods  have  been  used  to  solve  the  inverse  problem.  These 
methods  may  be  grouped  into  “direct”  and  “indirect”  methods  (Neuman,  1973: 1006).  The 
“direct”  methods  assume  the  model  parameters  are  dependent  variables  of  the  flow 
equation  (see  equation  3-13).  The  flow  equation  is  a  partial  differential  equation  of  the 
model  parameters.  Thus,  the  model  parameters  are  found  directly  by  solving  a  formal 
inverse  boundary  problem  (Carrera,  1987:559;  Yeh,  1986:96).  The  “indirect”  methods 
seek  to  improve  on  the  output  error  by  iteratively  adjusting  the  existing  parameter 
estimates  until  the  model  produces  results  “close”  to  the  real  system  (Neuman, 

1973:1006). 
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Direct  Methods.  To  apply  the  direct  approach,  the  hydrostatic  head  information 
must  be  known  over  the  entire  flow  region  and  the  measurement  errors  must  be  small.  In 
reality,  information  is  known  only  at  sparsely  distributed  wells.  To  properly  formulate  the 
inverse  problem,  additional  data  must  be  interpolated  to  fill  in  the  missing  nodal  points. 
With  the  aid  of  boundary  conditions  and  flow  data,  a  direct  solution  of  the  unknown 
parameters  in  the  original  governing  equation  may  be  made.  However,  when  the 
measured  and  interpolated  hydrostatic  head  values  are  entered  into  the  governing 
equation,  an  error  term  will  be  generated.  This  term  is  known  as  equation  error,  and  is 
minimized  by  the  proper  choice  of  parameters  (Yeh,  1986:96). 

Yeh  (1986:96)  and  Carrera  (1987:560-561)  list  a  number  of  investigators  that  use 
the  direct  method  to  solve  the  inverse  problem.  The  following  listing  presents  the  various 
techniques,  along  with  the  investigator: 


Energy  dissipation  (Nelson,  1968) 

Linear  programming  (Deininger,  1969;  Vemuri  and  Karplus,  1969; 
Kleinecke,  1971) 

Flatness  criterion  (Emsellem  and  de  Marsily,  1971) 
Multiple-objective  decision  process  (Neuman,  1973) 

Galerkin  method  (Frind  and  Pinder,  1973) 

Algebraic  approach  (Sagar  and  others,  1975) 

Inductive  method  (Nutbrown,  1975) 

Linear  and  quadratic  programming  (Hefez  and  others,  1975) 
Minimization  of  a  quadratic  objective  function 
with  penalty  function  (Navarro,  1977) 

Matrix  inversion  with  kriging  (Yeh  and  others,  1983) 


The  above  methods  are  not  based  on  a  statistical  framework.  The  variability  of  the 
parameters  and  the  error  terms  are  not  accounted  for  in  a  statistical  sense.  This  lack  of  a 
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statistical  foundation  may  be  the  reason  why  hydrologists  generally  do  not  use  the  above 
techniques  to  solve  the  inverse  problem  (Carrera,  1987:567).  Most  current  approaches  to 
solve  the  inverse  problem  are  statistically  based  and  can  be  categorized  as  “indirect” 
methods. 

Indirect  Methods.  The  indirect  methods  are  all  based  on  the  idea  of  minimizing 
the  difference  between  the  model’s  calculated  response  (typically,  hydraulic  head)  and  the 
actually  observed  value.  This  approach  is  well  suited  to  the  real-world  situation  where 
observation  wells  are  limited  in  number.  Because  the  relationship  between  hydraulic  head 
and  model  parameters  is  nonlinear,  indirect  methods  usually  require  an  iterative  approach 
to  solve  the  inverse  problem. 

Typically,  an  initial  set  of  parameters  is  chosen  for  the  groundwater  simulation 
model.  The  model’s  output  is  compared  to  the  field  observations  and  the  difference 
(error)  is  computed.  If  the  error  is  too  large,  one  or  more  parameters  are  adjusted,  and  the 
model  is  re-run.  If  this  error  is  within  a  pre-defined  performance  criterion,  the  latest 
parameters  are  adopted  for  the  model  (Neuman,  1973: 1006).  There  are:  1)  a  wide  variety 
of  ways  to  compute  the  error  between  the  modeled  and  measured  observations  and  2) 
numerous  methods  to  adjust  the  parameters  after  each  iteration.  These  characteristics  of 
the  indirect  approach  have  lead  to  numerous  publications  on  how  to  implement  the 
indirect  approach  to  solve  the  inverse  problem. 

A  manual  trial-and-error  method  is  most  often  used  to  arrive  at  estimates  of  the 
model  parameters.  This  method  is  labor  intensive  (therefore  expensive),  frustrating 
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(therefore  often  left  incomplete),  and  subjective  (therefore  biased)  (Carrera  and  Neuman, 
1986a:200). 


Yeh  (1986:96-98)  lists  a  number  of  investigators  that  used  indirect  methods  for 
parameter  estimation.  The  following  review  lists  the  various  techniques,  along  with  the 
investigator: 


Quasilinearization  and  control-oriented  techniques  (Bellman  and  Kabala,  1965; 

Yeh  and  Tauxe,  1971;  DiStefano  and  Rath,  1975) 

Minimax  and  linear  programming  (Yeh  and  Becker,  1973) 

Maximum  principle  (Lin  and  Yeh,  1974) 

Optimal  control  and  gradient  procedure  (Vermuri  and  Karplus,  1969) 

Optimal  control  using  steepest  descent  and  conjugate  gradient  method  (Chen  and 
others,  1974;  Chavent  and  others,  1975) 

Kalman  filtering  techniques  (McLaughlin,  1975;  Wilson  and  others,  1978) 
Maximum  likelihood  estimation  and  kriging  (Kitanidis  and  Vomvoris,  1983) 


Additional  inverse  problem  solution  methods  have  been  presented  since  the 
publication  of  Yeh  (1986).  The  papers  by  Carrera  and  Neuman  (1986a;  1986b;  1986c) 
advocate  using  a  maximum  likelihood  method  that  reaches  its  solution  by  means  of  a 
penalty  criterion  based  on  prior  estimates  of  the  model  parameters.  Xiang  and  others 
(1992)  use  an  LI  norm  to  measure  the  error  between  the  fitted  and  actual  response  in 
estimating  model  parameters.  Doughty  and  others  (1994)  base  their  method  on  iterated 
function  systems  using  fractal  theory. 

All  of  these  indirect  methods  primarily  deal  with  solving  the  inverse  problem  by 
calibrating  the  model  for  either  observed  head  or  contaminant  concentration  values. 
Usually,  the  flow  model  is  calibrated  first,  and  then  it  is  adjusted  to  properly  model  the 
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concentration  behavior.  Recently,  the  solution  to  the  coupled  flow-concentration  inverse 
problem  has  been  presented  in  the  literature  (Sun  and  Yeh,  1990;  Wagner,  1992). 


Coupled  Flow-Contaminant  Modeling 

The  common  method  of  estimating  the  hydrogeologic  parameters  independent  of 
the  characterization  of  the  contaminant  sources  requires  extra  iterations  to  eventually 
calibrate  the  model  for  both  flow  and  contaminants.  To  solve  for  the  hydrogeologic 
parameters,  the  contaminant  sources  must  be  assumed;  moreover,  to  estimate  the 
contaminant  sources,  the  hydrogeologic  parameters  must  be  assumed.  This  situation  can 
lead  to  numerous  trial-and-error  runs  until  both  objectives  are  calibrated. 

Wagner  (1992)  presents  a  method  that  solves  the  inverse  model  as  a  non-linear 
maximum  likelihood  problem.  The  method  simultaneously  estimates  the  hydrogeologic 
parameters  along  with  contaminant  source  characterization.  Due  to  the  statistical 
framework  of  this  maximum  likelihood  model,  statistical  estimates  of  the  uncertainties  of 
the  parameters  may  also  be  made. 

The  coupled  problem  may  also  be  applied  to  the  freshwater-saltwater  flow 
problem,  and  the  oil-water  two-phase  flow  problem  (Sun  and  Yeh,  1990:2508).  These 
problems  are  all  concerned  with  simultaneously  estimating  the  groundwater  model 
parameters  that  jointly  calibrate  the  model  for  two  types  of  response. 
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Summary 


The  published  approaches  to  calibrating  groundwater  flow  models  present 
numerous  alternatives  to  solving  an  ill-posed  inverse  problem.  Over  the  last  thirty  years, 
the  efficiency  and  precision  of  determining  the  model’s  parameters  have  improved.  A 
further  step  toward  efficient  calibration  of  groundwater  simulation  models  is  the 
simultaneous  estimation  of  the  parameters  for  both  groundwater  flow  and  contaminant 
transport.  However,  the  groundwater  model  calibration  literature  does  not  contain  an 
RSM  approach  to  solving  the  inverse  problem.  The  next  chapter  presents  how  the 
techniques  of  RSM  may  be  effectively  applied  to  groundwater  model  calibration. 
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III.  Groundwater  Model  Calibration 


Introduction 

This  chapter  reviews  the  fundamental  principles  involved  in  groundwater  model 
calibration.  It  first  presents  the  geologic  framework  and  hydrogeological  principles  of 
groundwater  modelling.  These  principles  are  then  used  to  develop  the  governing 
differential  equation  found  in  most  numerical  groundwater  flow  models.  Finally,  the 
problem  of  calibrating  the  groundwater  flow  model  and  how  response  surface  methods 
may  be  applied  to  the  calibration  effort  are  discussed. 


Groundwater  Hydrogeology 

Geological  Terms.  The  mathematical  model  that  represents  the  physical 
processes  of  a  groundwater  system  is  based  upon  the  geological  properties  of  the  rocks 
and  sediments  in  the  modeled  area.  A  basic  understanding  of  these  properties  and  the 
associated  terminology  is  necessary  to  develop  the  governing  equation  of  the  groundwater 
flow  model. 

The  aquifer  is  the  geological  unit  of  porous  material  capable  of  storing  and 
transmitting  appreciable  quantities  of  water  to  wells  (Anderson  and  Woessner,  1992: 12). 
The  aquifer  is  often  composed  of  several  different  rock  and  sediment  types  which  vary 
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spatially.  This  heterogeneous  nature  of  aquifers  is  often  modeled  by  dividing  the  aquifer 
into  zones  which  posses  similar  hydraulic  properties. 

The  aquifer  itself  is  composed  of  solid  matter  called  the  solid  matrix  and  void 
spaces  between  individual  grains  or  in  cracks  (fractures),  which  are  collectively  referred  to 
as  the  pore  space.  Fractures  are  the  separated  cracks  that  exist  in  many  massive  rock 
formations.  When  the  inter-granular  pore  space  has  been  reduced  or  eliminated  by 
minerals  and  cements,  fractures  may  provide  the  dominant  storage  space  in  an  aquifer.  In 
fractured  aquifers,  the  frequency  and  direction  of  the  fractures  will  control  the  amount  and 
ease  of  movement  of  the  water  from  one  area  to  another. 

Hydraulic  Properties.  There  are  several  properties  of  aquifers  that  characterize 
groundwater  systems.  The  properties  of  most  importance  to  groundwater  modeling  are 
porosity,  hydraulic  conductivity,  and  permeability. 

Porosity  is  the  ratio  of  the  volume  of  pore  space  to  the  total  volume  of  the  aquifer 
system.  This  ratio  is  a  dimensionless  number  less  than  1,  but  most  rocks  have  a  porosity 
below  0.30.  Fractured  rocks  exhibit  a  wide  range  of  porosity  from  0.01  to  0. 10  .  In  the 
saturated  zone  of  the  aquifer,  the  water  resides  in  this  pore  space,  and  the  amount  of  water 
contained  per  unit  volume  is  directly  measured  by  the  porosity.  For  groundwater  flow 
systems,  only  the  interconnected  pores  that  permit  the  flow  of  water  are  considered  (Smith 
and  Wheatcraft,  1992:6.9).  This  effective  porosity  determines  the  amount  of  water  per 
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unit  volume  that  is  available  for  flow  and  is  defined  as 


Effective  Porosity  = 


Interconnected  Void  Volume 
Total  Volume 


(3-1) 


Hydraulic  conductivity  measures  the  ability  of  a  fluid  to  move  through  a  porous 
medium.  The  hydraulic  conductivity  ( K)  is  a  function  of  both  the  medium  and  the  fluid 
and  is  defined  as 


kP-P-S- 

P 


(3-2) 


where  k  (m2)  is  the  permeability,  a  property  of  the  medium  only  that  describes  how  well 
the  medium  transmits  water;  p  (kg /  m3)  is  the  fluid  density;  g  (m/min2)  is  the  acceleration 
due  to  gravity;  and  p  (kg/m-min)  is  the  dynamic  viscosity  of  the  fluid,  giving  K  units  of 
m/min  (Smith  and  Wheatcraft,  1992:6.9). 


Groundwater  Flow 

Hydraulic  Head.  Under  the  influence  of  gravity,  groundwater  flows  from  regions 
of  high  hydraulic  head  to  regions  with  lower  hydraulic  head.  Hydraulic  head  (h)  is 
measured  in  meters  and  is  defined  as 


h-  z  + 


P 

P  g 


(3-3) 
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where  p  (g  kg/m2)  is  the  pressure  of  the  fluid  with  constant  density  p,  and  g  is  the 
acceleration  due  to  gravity.  2  is  the  elevation  of  the  measure  point  above  datum  (a  level 
reference  height),  and  is  known  as  the  elevation  head. 

The  water  table  is  defined  as  the  surface  on  which  the  pressure  head  (hp)  is  equal 
to  zero.  The  pressure  head  is  expressed  in  terms  of  values  above  atmospheric  pressure  as 
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The  relationship  among  hydraulic  head,  pressure  head  ,  and  elevation  head  are  illustrated 
in  Figure  3.1  (Smith  and  Wheatcraft,  1992:6.7). 

Figure  3.1  Relationship  Among  Heads 
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Groundwater  Flow  Governing  Equation.  The  basic  principles  of  groundwater 
flow  were  investigated  and  reported  by  Henri  Darcy  in  1 856.  His  studies  established  that 
the  flow  rate  through  porous  media  is  proportional  to  the  hydraulic  head  loss  and  inversely 
proportional  to  the  length  of  the  flow  path.  Darcy’s  law  in  three  dimensions  is 
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or 

q  =  -KVh  (3-5) 

where  q  is  the  specific  discharge  rate  (m/min),  K  is  the  hydraulic  conductivity,  and  h  is  the 
hydraulic  head. 

To  derive  the  general  form  of  the  governing  equation  of  groundwater  flow,  first 
consider  a  representative  elementary  volume  (REV)  of  the  aquifer,  as  shown  in  Figure  3.2. 


Figure  3.2  Representative  Elementary  Volume  (Anderson  and  Woessner,  1992:16) 


This  cube  represents  a  portion  of  the  aquifer  large  enough  to  represent  the  properties  of 
the  porous  medium,  yet  small  enough  that  the  change  in  head  within  this  cube  is  relatively 
small.  The  volume  of  this  REV  is  AxAyAz  .  The  flow  of  water  through  this  volume  is 
expressed  in  terms  of  the  specific  discharge  rate  q,  which  may  be  expressed  as 


q  =  q£x  +  <?A  +  (3  -6) 

in  terms  of  the  orthogonal  unit  vectors  ix,  iy,  and  iz  along  the  x ,  y,  and  z  axes,  respectively. 
By  conservation  of  mass1,  the  rate  of  change  in  storage  is  defined  as 

Outflow  rate  -  Inflow  rate  =  Rate  Of  Change  in  Storage  ( ROCS)  (3-7) 

Considering  only  the  flow  along  the  y  axis,  the  change  in  flow  rate  through  the  REV  is 

d  qv 

— - ^-(AxAyAz)  (3-8) 

d  y 


Similar  expressions  may  be  written  for  the  x  and  z  axes.  The  sum  of  the  change  in  flow  for 
all  three  axes,  minus  a  term  for  the  recharge  R,  is  equal  to  the  ROCS: 


rs  qx  ,  g  qy 

v  d  x  d  y 


C-^-R 

d  z 


AxAyAz  =  ROCS 


(3-9) 


The  ROCS  is  the  ratio  of  change  in  volume  per  unit  of  time.  The  change  of 
volume  is  described  by  specific  storage  Ss,  defined  as  the  volume  of  water  released  from 
storage  per  unit  change  in  head  h ,  per  unit  volume  of  aquifer: 


1  With  the  assumption  of  constant  density,  conservation  of  volume  is  also  maintained. 
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(3-10) 


S=- 


AV 


AhAxAyAz 


Therefore,  the  rate  of  change  in  storage  in  the  REV  is 


ROCS  =  ~  =  ~Ss—  Ax  A yAz 
At  ‘A  t  y 


(3-11) 


By  combining  equations  3-9  with  3-11,  dividing  through  by  AxAyAz ,  and  taking 
the  limit  as  At  -»  0 ,  the  final  theoretical  form  of  the  water  balance  equation  is 


^L+13l+£JL  =  -S,^  +  r 

d  x  8  y  8  z  8 1 


(3-12) 


Because  the  discharge  rate  q  cannot  be  measured  directly,  this  equation  cannot  be  applied 
to  field  observations.  However,  Darcy’s  law  (equation  3-5)  defines  the  relation  between  q 
and  /z;  and  head  is  a  parameter  that  may  be  directly  measured  in  the  field.  This  equation 
may  be  further  simplified  for  steady-state  simulations.  Under  steady-state  conditions,  the 

8  h 

heads  will  not  change  with  time,  thus -  equals  zero.  If  we  additionally  assume  that  no 

8  t 

recharge  occurs,  then,  when  the  equations  of  3-5  are  substituted  into  equation  3-12,  the 
governing  equation  of  groundwater  flow  under  steady-state  conditions  is 
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(3-13) 


Some  form  of  equation  3-13  is  used  in  all  steady-state,  numerical,  groundwater  flow 
models  with  no  recharge  (Anderson  and  Woessner,  1992:13-18). 
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Numerical  Groundwater  Modeling. 


The  numerical  model  SUTRA  (Saturated- Un saturated  TRAnsport)  was  used  to 
simulate  groundwater  flow  for  this  study.  SUTRA  was  developed  by  the  United  States 
Geological  Survey  for  computer  modeling  of  groundwater  systems.  The  version  used  in 
this  study,  V-0690-2D,  was  released  in  June  of  1990.  The  program  employs  a  two- 
dimensional  hybrid  finite-element  and  integrated  finite-difference  method  to  approximate 
the  governing  differential  equation  presented  as  equation  3-13  (Voss,  1984:3). 

The  numerical  method  used  by  SUTRA  divides  the  cross-section  of  the  aquifer 
into  a  single  layer  of  contiguous  blocks,  called  finite  elements.  Flow  parameters  and 
variables  which  vary  from  point  to  point  within  the  aquifer  are  approximated  at  these 
elements.  The  comers  of  the  elements  are  termed  nodes,  and  regions  that  are  centered  at 
nodes  are  called  cells  (figure  3.3).  This  grid  of  interconnected  nodes  may  be  constructed 
on  a  cross-sectional  view  of  the  aquifer  being  modeled.  The  spacing  and  placement  of  the 
nodes  will  depend  on  the  actual  geology  of  the  aquifer  and  will  need  to  be  adjusted  by  the 
analyst  to  accurately  reflect  the  hydrogeological  changes  of  the  aquifer. 

Figure  3.3  Graphical  Representation  of  Elements,  Nodes,  and  Cells 
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Each  execution  of  SUTRA  requires  an  input  parameter  file  to  specify  the  hydraulic 
conductivities  and  porosity  for  each  element  of  the  grid.  Once  this  parameter  file  is 
prepared,  SUTRA  solves  equation  3-13  and  provides  the  steady-state  fluid  pressures  at 
every  node  of  the  grid  as  its  primary  output.  When  SUTRA  is  run  in  steady-state  mode, 
the  fluid  pressures  do  not  vary  with  time.  For  convenience,  these  output  pressures  can  be 
converted  to  hydraulic  head  ( h )  using  equation  3-3.  Appendixes  B  and  C  contain  the 
FORTRAN  code  and  VAX  command  files  that  were  developed  to  streamline  the 
execution  of  SUTRA  for  this  project.  These  files  allowed  for  rapid,  interactive  data  entry 
and  automatic  reporting  of  results. 


Calibration 

As  mentioned  above,  SUTRA  requires  several  parameters  for  each  of  the  elements 
of  the  grid  to  solve  the  finite-element  model  of  the  aquifer.  The  process  of  adjusting  the 
model’s  parameters  (e.g.  hydraulic  conductivity  and  effective  porosity)  to  obtain  a 
reasonable  but  nonunique  match  between  the  observed  site-specific  data  (calibration 
target)  and  the  model’s  estimated  values  is  known  as  model  calibration  (Walton, 

1992:35).  The  difference  between  an  observed  calibration  target  and  the  model’s  estimate 
of  that  target  is  the  residual  or  error  in  the  estimate. 

The  calibration  target  dataset  for  this  study  was  composed  of  the  hydraulic  head 
values  ( hm )  observed  (measured)  within  an  aquifer.  By  computing  the  differences  between 
the  heads  in  the  calibration  target  ( hm )  with  the  heads  computed  from  the  finite-element 
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model  of  SUTRA  (h,\  the  accuracy  of  the  model  may  be  assessed.  For  this  study,  four 
residual  statistics  were  computed  for  each  execution  of  the  SUTRA  model,  where  the 
number  n  of  elements  in  the  model  was  524.  The  equations  for  calculating  the  Sum  of 
Squared  Error  (SSE),  Root  Mean  Squared  Error  (RMS),  Mean  Absolute  Error  (MAE), 
and  Mean  Error  (ME)  are  respectively  (Anderson  and  Woessner,  1992:238-241) 


n 


SSE  =  £  (K-h,)] 

2=1 

(3-14) 

r  i.  i0-5 

2=1  _ 

(3-15) 

2=1 

(3-16) 

(3-17) 

2=1 


During  the  calibration  process,  the  input  parameters  to  the  SUTRA  model  are 
adjusted  until  a  suitably  small  residual  statistic  is  computed  using  one  of  the  above 
equations.  As  discussed  in  Chapter  II,  numerous  methods  have  been  proposed  to 
analytically  determine  the  optimal  choice  of  input  parameters  that  will  provide  a  calibrated 
model  with  a  small  error.  Response  Surface  Methodology  (RSM)  provides  a  means  for 
systematically  altering  the  input  parameters  to  the  SUTRA  model  to  achieve  a  suitably 
small  residual. 
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Response  Surface  Methodology 


Response  surface  methodology  comprises  a  group  of  statistical  techniques  for 
empirical  model  building  and  exploitation.  One  primary  goal  of  RSM  is  to  find  the  best 
value  of  the  estimated  response.  The  response  chosen  for  this  study  was  SSE  (equation  3- 
14)  because  SSE  was  believed  to  be  amenable  to  local  approximation  by  the  low-order 
polynomials  typically  used  in  RSM  studies.  The  goal  would  be  to  find  the  set  of  parameter 
values  which  minimize  the  SSE,  thus  yielding  calculated  hydraulic  heads  from  the  finite- 
element  groundwater  model  that  closely  match  the  calibration  data  set. 

Using  RSM  to  calibrate  a  groundwater  model  involves  three  basic  phases  of 
investigation.  The  first  phase  is  a  screening  process  that  determines  which  parameters 
significantly  influence  the  value  of  the  estimated  SSE  response.  The  second  phase  uses 
local  first-order  approximations  of  the  SSE  response  surface  to  determine  parameter 
values  associated  with  a  smaller  SSE.  Finally,  the  third  phase  employs  second-order 
approximations  of  the  SSE  response  surface  to  identify  locally  optimal  parameter  settings 
that  minimize  the  SSE  for  the  calibrated  model. 

A  key  feature  of  the  RSM  technique  is  the  empirical  testing  of  the  response 
surface.  At  no  time  is  an  attempt  made  to  determine  the  shape  of  the  entire  response 
surface;  rather,  through  first-  and  second-order  models,  pieces  of  the  response  surface  are 
approximated  and  actual  tests  made  in  the  direction  of  suggested  improvement.  This 
empirical  testing  leads  to  an  understanding  of  the  underlying  SSE  response  without  the 
need  to  precisely  define  it  (Box  and  Draper,  1987:1-19). 
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I  Parameter  Screening  Phase.  The  first  step  in  calibrating  a  groundwater  model 
is  to  determine  which  parameters  need  to  be  adjusted.  Parameters  such  as  porosity  and 
hydraulic  conductivity  may  be  adjusted  for  as  many  hydrogeologic  zones  as  are  deemed 
necessary  to  account  for  the  presumed  differences  in  the  aquifer.  This  approach  may  lead 
to  a  dozen  or  more  parameters  that  need  to  be  set  in  the  groundwater  model.  Fortunately, 
only  relatively  few  of  these  parameters  are  likely  to  significantly  influence  the  behavior  of 
the  groundwater  model.  The  purpose  of  the  screening  experiment  is  to  determine  which 
parameters  have  a  significant  effect  on  the  groundwater  model’s  estimate  of  the  hydraulic 
heads.  Thus,  the  screening  phase  serves  the  same  purpose  as  sensitivity  analysis  in 
traditional  groundwater  modeling  studies. 

The  parameter  screening  phase  usually  relies  on  a  two-level  experimental  design. 

In  a  two-level  design,  the  parameters  are  set  to  high  and  low  levels  of  their  expected 
range,  based  on  the  observed  local  hydrogeological  conditions.  For  convenience  and 
computational  efficiency,  these  parameter  values  are  often  linearly  transformed  into 
coded  or  standardized  variables  xk,  such  that 

—  (3-18) 

where  E,k0  is  the  center  of  the  region  of  interest  and  Sk  is  the  half  range  of  the  region.  A 
two-level  factorial  (2k)  design  where  each  of  the  k  variables  occurs  at  just  two  levels  (±1  in 
coded  space)  may  be  used  to  estimate  the  k  main  effects  in  2k  runs.  As  an  alternative 
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especially  suited  for  parameter  screening,  Plackett-Burman  designs  allow  for  estimation  of 
the  k  main  effects  in  only  k  +  1  runs  (Box  and  Draper,  1987: 162). 

Once  the  experimental  runs  have  been  performed,  the  SSE  responses  are  fit  to  a 
first-degree  polynomial  model  in  k  coded  variables  (Cornell,  1990: 13) 

K=  Po+  P\Xu\  +  PlXu2  +•  •  ‘+PkXuk  +  £u  (3-19) 

The  coefficients  are  proportional  to  the  effect 'the  Ath  variable  has  on  the  output  of  the 
groundwater  model  (Effect^  =  2 /3k).  The  presumption  of  the  screening  phase  is  that  only 
those  parameters  with  significant  first-order  effects  need  to  be  considered  in  subsequent 
experiments.  Normally,  ANOVA  methods  would  indicate  which  coefficients  are 
significant,  but  the  relatively  few  runs  in  a  Plackett-Burman  design  preclude  this  approach. 
In  such  instances,  the  significant  effects  can  be  identified  through  use  of  a  normal 
probability  plot. 

In  a  properly  fit  first-order  linear  model,  the  residuals  are  approximately  normally 
distributed  with  equal  variance.  If  there  are  no  significant  effects,  each  of  the  computed 
effects  will  represent  an  observation  from  a  common  normal  error  distribution.  A  plot  on 
normal  probability  paper  of  the  empirical  cumulative  distribution  function  of  these 
ordered,  insignificant  effects  should  roughly  form  a  straight  line.  However,  if  some  of  the 
plotted  points  differ  noticeably  from  a  straight  line  either  at  the  top  right  or  bottom  left, 
then  the  corresponding  effects  are  presumed  to  be  significant. 

Standard  graph  paper  may  be  used  to  identify  the  significant  effects  using  the 
following  procedure: 
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1)  Order  the  effects  from  smallest  to  largest  and  scale  the  x-axis  accordingly. 


2)  Plot  the  quantity  ^[(i-O.Syk],  where  (jf1  (p)  is  the  inverse  cumulative 
distribution  function  of  the  standard  normal  distribution  and  i  represents  the  rank  of  the 
effect.  This  function  may  be  approximated  by 

f1  (p)  =  [p0  1349  -  (l-p)01349]/0.1975  (3-20) 

From  this  plot,  only  those  parameters  which  significantly  contribute  to  the  output 
of  the  groundwater  model  are  selected  for  further  analysis.  In  the  next  two  phases,  the 
settings  of  these  parameters  are  adjusted  to  reduce  the  value  of  the  SSE  response  to  within 
calibration  tolerances. 

II  First-Order  Design  Phase.  Response  surface  methodology  investigates  the 
nature  of  the  SSE  response  surface  by  conducting  designed  experiments  centered  at 
particular  design  points.  A  two-level,  full-factorial  (2*)  design,  where  each  of  the  k 
variables  occurs  at  just  two  levels  (±1  in  coded  space),  can  be  used  to  empirically  evaluate 
the  SSE  response  surface.  This  design  requires  2k  runs.  Using  the  responses  from  the 
actual  process  (the  SSE  responses  from  individual  SUTRA  runs),  an  estimate  of  the 
gradient  of  the  SSE  response  surface  in  the  design  region  may  be  made.  The  estimate  of 
the  gradient  is  determined  by  fitting  a  first-order  model  like  equation  3-19  to  estimate  the 
SSE  responses  actually  observed  as  a  function  of  the  coded  input  values.  This  first-order 
model  of  the  response  surface  is  analogous  to  using  a  truncated  Taylor’s  series  to 
approximate  a  function. 
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Once  the  gradient  direction  is  determined  (from  the  coefficients  of  the  first-order 
model),  further  experiments  are  conducted  in  the  direction  of  steepest  descent.  The 
experiments  are  conducted  away  from  the  center  of  the  starting  design  region  until  the 
response  begins  to  increase  (in  groundwater  model  calibration,  the  goal  is  to  locate  a 
minimum  in  the  SSE  response).  Once  a  minimum  along  a  gradient  is  reached,  a  new 
design  may  be  established,  a  new  gradient  direction  computed,  and  the  gradient  search 
along  the  steepest  descent  may  be  continued. 

The  first-order  designs  may  be  iteratively  executed  as  long  as  measures  to  asses  the 
adequacy  of  fit  indicate  a  first-order  design  is  still  appropriate.  A  first-order  model  should 
have  good  explanatory  power.  This  explanatory  power  is  demonstrated  by  relatively  high 
R  Square  and  F  values  of  the  regression  used  to  determine  the  coefficients  of  the  first- 
order  model.  Also,  a  reasonable  amount  of  improvement  in  the  response  should  be 
realized  over  the  previous  design  point.  Finally,  a  test  for  curvature  should  indicate  there 
is  no  significant  curvature  in  the  fitted  first-order  model.  This  test  evaluates  whether  a  run 
performed  at  the  design  center  falls  on  a  hyperplane  running  through  the  2k  factorial  points 
or  if  a  significant  amount  of  curvature  is  needed  to  fit  the  center  point.  Evidence  of 
curvature  implies  significant  higher-order  terms  have  been  omitted. 

Single  Degree  of  Freedom  Test  for  Curvature  In  order  to  statistically  evaluate  the 
adequacy  of  the  first-order  model,  the  center  point  run  can  be  used  in  a  Single  Degree  of 
Freedom  Test  for  Curvature.  A  Sum  of  Squares  of  Pure  Quadratic  terms  ( SSPO )  is  given 
by  the  following  “curvature  contrast”: 
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where  n0  is  the  number  of  center  point  runs  (which  equals  1  in  the  case  of  a  deterministic 
model  like  SUTRA),  nf  is  the  number  of  factorial  point  runs  (which  equals  2k  for  a  full 

factorial  design),  y  is  the  average  of  the  SSE  responses  for  the  factorial  points,  and 

is  the  SSE  response  observed  at  the  center  point  run.  The  Sum  of  Squares  Residual 

( SSR )  presented  in  an  ANOVA  table  is  partitioned  into  Quadratic  Terms  and  Other  Terms 
as  follows: 

Quad.  Terms  =  SSPQ  (3-22) 

Other  Terms  =  SSR  -  SSPQ  (3-23) 

These  two  terms  are  divided  by  their  respective  degrees  of  freedom  forming  Mean  Square 
of  Quadratic  (MSQ)  and  Other  Terms  ( MSO ).  The  ratio  of  MSQ/MSO  is  an  F  statistic 
which  may  be  compared  to  an  F-distribution  with  1  and  2*-k-l  degrees  of  freedom,  for  a 
full  factorial  design,  corrected  for  the  mean.  A  high  F-value  or  the  correspondingly  low  p- 
value  implies  significant  quadratic  effects.  Thus,  when  this  test  fails,  a  first-order  model 
(equation  3-19)  should  probably  not  be  used  to  approximate  the  response  surface  at  that 
design  point.  In  general,  if  several  measures  indicate  the  first-order  model  is  not  adequate, 
the  first-order  phase  should  be  terminated  and  the  second-order  design  phase  should  be 
pursued. 


HI  Second-Order  Design  Phase.  The  primary  reason  the  first-order  strategy  had 
to  be  abandoned  was  because  the  first-order  polynomial,  equation  3-19,  was  not  able  to 
adequately  fit  the  observed  responses.  In  effect,  the  truncated  Taylor’s  series  had  been 
truncated  too  far.  In  the  second-order  design  phase,  the  following  second-degree 
polynomial  model  is  used  (Box  and  Draper,  1987:376): 

Y  —  P 0  +  A*,  ^  ^  P  kXk  PnXU^  *" PkkXk  +  P\2X\X2  4  ^  P  k-\,kX  k-\X  k  (3-24) 

The  two-level,  full-factorial  (2*)  design  used  in  the  first-order  phase  is  not  able  to  provide 
enough  information  to  estimate  the  quadratic  terms  of  equation  3-24.  To  build  a  quadratic 
model,  a  central  composite  design  (CCD)  may  be  used.  A  central  composite  design 
consists  of  the  following  three  parts 

i)  the  2k  vertices  (for  a  full  factorial  design)  of  a  ^-dimensional  cube  where  the 
factor  levels  are  coded  so  that  the  design  center  is  at  (0,  0, .  . .  ,0).  The  values 
of  the  coded  variables  in  this  factorial  portion  are  (±1,  ±1, .  .  . ,  ±1); 

ii)  the  center  point,  no=(0,  0, . . . ,  0)  and; 

iii)  the  2k  vertices  (±a,  0,  0, ... ,  0),  (0,  ±a,  0, .  .  . ,  0), .  .  .  (0,  0, .  .  . ,  0,  ±a)  of  a 
^-dimensional  “star”  (Cornell,  1990:52).  Figure  3.4  shows  a  central  composite 
design  for  three  variables  ( k  =  3). 
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Figure  3.4  CCD  for  k  =  3  (Cornell,  1990:53) 


A  desirable,  but  not  necessary  feature  often  incorporated  into  a  CCD  is 
rotatability.  In  a  rotatable  design,  the  accuracy  of  prediction  of  the  response  depends  only 
on  its  distance  from  the  center  of  the  design,  producing  a  variance  function  that  is  nearly 
spherical.  Rotatability  is  achieved  by  assigning  appropriate  a  values  to  the  “star”  points  of 
the  CCD.  A  general  formula  used  to  achieve  rotatability  in  a  full  factorial  design  is: 

a  =  (: 2kf  (3-25) 

where  k  is  the  number  of  factors  (Cornell,  1990:52).  For  a  more  general  formula,  see  Box 
and  Draper  (1987:488). 

One  of  the  advantages  of  the  central  composite  design  is  that  the  factorial  points  in 
(i)  above  are  the  same  points  used  in  the  last  two-level,  full-factorial  (2h)  design  from  the 
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first-order  phase.  In  most  cases,  the  center  point  response  would  already  have  been 
determined  from  either  the  single  degree  of  freedom  test  for  curvature  or  the  steepest 
descent  gradient  search  conducted  from  the  center  of  the  last  first-order  design.  Thus,  to 
begin  the  second-order  design  phase,  only  the  star  points  are  run  to  obtain  the  additional 
responses  needed  to  estimate  the  coefficients  of  equation  3-24. 

Once  this  second-order  model  has  been  fit  to  the  data,  the  method  of  steepest 
descent  applied  to  second-order  surfaces  needs  to  be  performed  to  identify  potential 
design  points  to  use  in  experimenting  beyond  the  design  region.  This  second-order 
surface  method  of  steepest  descent  is  known  as  ridge  analysis. 

Ridge  Analysis.  The  method  of  ridge  analysis  locates  the  minimum  estimated 
response  on  the  fitted  model,  equation  3-24,  that  is  a  distance  R  from  the  design  center. 
Ridge  analysis  may  be  summarized  as  the  following  constrained  optimization  problem: 

Minimize  Y  =  fio  +  x'fi+x'fix 
Subject  to  xf  -  R2  =  0 

i 

This  optimization  problem  can  be  solved  using  the  SAS  statistical  package  to 
automatically  exploit  a  ridge  system.  By  conducting  a  ridge  analysis  for  each  of  several 
distances  R  from  the  design  center,  SAS  provides  a  set  of  promising  points  to  use  in 
establishing  the  location  of  improved  responses.  Ridge  analysis  does  not  conduct  a  search 
along  the  specified  path;  it  only  identifies  potential  design  points  to  use  in  experimenting 
beyond  the  design  region  (Box  and  Draper,  1987:375-380).  By  actually  conducting 
experiments  along  the  identified  path,  the  responses  will  typically  improve  (decrease)  and 
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then  become  worse  (increase)  as  the  second-order  model  diverges  from  the  true 
underlying  response  surface.  If  the  lowest  response  found  provides  a  groundwater  model 
which  is  calibrated  within  tolerances,  the  hydrologist  is  finished.  If  further  improvement  is 
sought,  another  CCD  may  be  established  at  the  lowest  response,  and  the  second-order 
design  phase  repeated. 

The  iterative  use  of  response  surface  methods  will  identify,  within  error  tolerances, 
a  local  minimum  of  the  response  surface.  The  parameters  of  the  groundwater  model  used 
at  this  minimum  have  calibrated  the  finite-element  model.  Although  several  first-  and 
second-order  experimental  designs  may  be  needed  to  locate  this  minimum,  a  suitably 
calibrated  model  could  possibly  be  achieved  before  this  local  minimum  is  found.  In  this 
case,  the  hydrologist  may  elect  to  discontinue  the  iterative  procedure  once  the  desired 
degree  of  tolerance  is  achieved. 
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IV.  Results  and  Analysis 


Introduction 

The  primary  objective  of  this  project  was  to  determine  if  the  techniques  of 
Response  Surface  Methodology  (RSM)  could  be  used  to  calibrate  a  groundwater  flow 
model.  First,  a  set  of  steady-state  hydraulic  head  values  was  established  to  validate  the 
proper  functioning  of  the  finite-element  code  being  used  (SUTRA,  Saturated- Unsaturated 
Z&4nsport).  These  steady-state  hydraulic  head  values  also  served  as  the  calibration  target 
data  set  to  compare  against  the  flow  model’s  output.  Next,  a  screening  experiment  was 
conducted  to  determine  which  parameters  had  a  significant  effect  on  the  flow  model’s 
outcome.  Finally,  the  RSM  approach  was  used  to  calibrate  this  two-dimensional,  steady- 
state  flow  model  by  adjusting  the  significant  parameters  until  a  local,  if  not  global, 
optimum  was  obtained. 

Calibration  Target  Data  Set  Preparation  /  Validation 

In  choosing  an  example  to  validate  the  RSM  approach,  a  previously  calibrated 
groundwater  model  was  used  as  a  calibration  target  rather  than  calibrating  a  new  model 
using  field  observations.  Using  a  previously  calibrated  model  also  had  the  advantage  of 
eliminating  uncertainties  due  to  field  measurement  errors.  This  approach  effectively 
provided  observational  head  data  at  each  computational  node  in  the  calibration  target  data 
set.  While  such  resolution  is  unrealistic  for  field  applications,  it  nonetheless  provided  an 
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excellent  platform  for  testing  the  RSM  procedure.  This  validation  approach  is  commonly 
seen  in  the  literature  (Xiang  and  others,  1993;  Willis  and  Yeh,  1987;  Carrera  and  Neuman, 
1986c). 

Smith  and  Ritzi  (1993)  used  a  two-dimensional,  finite-element  model  in  cross- 
section  mode  to  simulate  groundwater  flow  and  nitrate  transport  on  the  Sycamore  Farm 
research  facility  of  Wright  State  University,  Ohio.  A  complete  hydrogeologic  description 
of  the  Sycamore  Farm  area  is  provided  in  Smith  (1991 :55-56).  The  Smith-Ritzi  model 
served  as  the  calibration  target  data  set  for  this  thesis  effort.  As  the  input  parameters  were 
adjusted  in  the  SUTRA  model  using  the  principles  of  RSM,  the  output  was  compared  to 
the  Smith-Ritzi  model  to  determine  how  closely  the  SUTRA  model  was  being  calibrated. 

The  Smith-Ritzi  model  contains  ten  hydraulic  conductivity  zones.  The  values  used 
for  these  ten  hydraulic  conductivity  zones  are  presented  in  table  4. 1 . 


Table  4.1  Hydraulic  Conductivity  Zones 


Hydraulic  Conductivity 
Zone 


Unit  1 
Unit  2 
Unit  3 
Unit  4 
Unit  5 
Unit  6 
Unit  7 
Unit  8 
Unit  9 
Unit  10 


Geological  Formation 


Dayton  Dolomite 
Upper  Brassfield 
Upper  Brassfield 
Upper  Brassfield 
Upper  Brassfield 
Gray  Clayey  Shale 
Middle  Brassfield 
Middle  Brassfield 
Lower  Brassfield 
Lower  Brassfield 


Hydraulic  Conductivity 
(m/min) 


2.0  x  10'3 
1.0  x  10'3 
1.0  x  10-4 
3.0  x  10'5 
2.0  x  10’3 
3.0  x  lO"4 
6.0  x  10-4 
6.0  x  10-6 
4.0  x  10-6 
1.0  x  10’5 
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The  hydraulic  conductivities  listed  above  and  a  single  effective  porosity  of  1 1% 
were  used  as  the  parameters  for  the  SUTRA  program  to  produce  a  calibration  target  set  of 
heads.  This  raw  output  contained  the  horizontal  (X)  and  vertical  (Y)  grid  coordinates  as 
well  as  the  hydraulic  pressure  for  each  of  the  524  nodes  in  the  finite-element  grid.  Figure 
4. 1  shows  the  finite-element  grid  used  by  SUTRA  to  compute  the  hydraulic  pressures  at 
the  524  nodes.  The  aquifer  being  modeled  consists  of  nine  stratigraphic  units,  divided  into 
ten  hydrologic  units,  each  with  its  own  hydraulic  conductivity.  The  spacing  and 
arrangement  of  these  nodes  was  determined  by  the  local  geology  of  the  groundwater 
system  and  reflects  the  hydraulic  conductivity  zonation  shown  in  figure  4.2. 

The  raw  output  from  the  SUTRA  run  was  not  in  units  that  were  directly  useful. 

The  vertical  coordinates  (Y)  were  given  in  meters  above  a  horizontal  datum  fixed  at 
274.35  meters  above  mean  sea  level  (AMSL).  These  vertical  coordinates  were 
transformed  to  heights  AMSL  (Yamsl)  by  simply  adding  274.35  m  to  all  524  vertical 
coordinates.  The  hydraulic  pressures  computed  by  SUTRA  at  each  node  were  transformed 
to  steady-state  hydraulic  heads  via  the  following  function: 

HeadAMSL  =  Pressure/9800  +  Y  amsl  (4- 1 ) 

where  the  pressure  is  given  in  units  of  kg/m- s2  and  heights  AMSL  in  terms  of  m. 
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Figure  4.1  Cross-Sectional  Finite-Element  Grid  (Smith,  1991:28) 
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After  completing  the  above  calculations,  a  calibration  target  data  set  consisting  of 
the  steady-state  hydraulic  heads  at  all  524  nodes  was  ready.  It  contained  the  X  and  Yamsl 
coordinates  in  meters  from  the  lower  left  (southwest)  comer  of  the  cross  section,  along 
with  their  accompanying  steady-state  hydraulic  head  values  (HeadAMSL),  also  measured  in 
meters  AMSL.  This  calibration  target  data  set  of  hydraulic  heads  was  in  perfect 
agreement  with  the  hydraulic  heads  presented  by  Smith  (1991 :35)  (see  figure  4.3).  By 
computing  a  set  of  heads  that  matched  the  heads  of  the  calibrated  Smith-Ritzi  model,  the 
calibration  target  set  was  validated  for  use  in  the  following  calibration  study.  All 
subsequent  outputs  from  the  SUTRA  model  were  compared  to  this  calibration  target  set 
to  determine  the  progress  of  the  calibration. 
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Screening  Experiments  /  Sensitivity  Analysis 

The  aquifer  being  modeled  contained  ten  hydraulic  conductivity  zones,  and  an 
overall  effective  porosity.  Thus,  eleven  possible  parameters  may  be  adjusted  to  try  and 
match  the  output  of  the  SUTRA  model  with  the  calibration  target  set  of  hydraulic  heads. 
Before  attempting  to  adjust  1 1  parameters  to  calibrate  the  model,  a  screening  experiment 
or  sensitivity  analysis  was  conducted  to  determine  which  parameters  significantly 
contributed  to  the  model’s  output. 

A  screening  experiment,  based  on  a  two-level  Plackett-Burman  design  (Box  and 
Draper,  1987:162),  was  established  to  evaluate  the  significance  of  the  11  parameters. 
Table  4.2  presents  the  settings  used  for  the  12  runs  of  this  screening  experiment.  The 
complete  results  are  listed  in  Appendix  A. 


Table  4.2  Plackett-Burman  Experimental  Values 


Run 

Unit  1 

Unit  2 

Unit  3 

Unit  4 

Unit  5 

Unit  6 

Unit  7 

Unit  8 

Unit  9 

Unit  10 

Porositv 

1 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

20 

2 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1 

3 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

20 

4 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

20 

5 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

20 

6 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1 

7 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-07 

1 

8 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1 

9 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-07 

20 

10 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

1.0E-02 

1 

11 

1.0E-07 

1.0E-02 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-02 

1.0E-02 

1.0E-02 

1.0E-07 

1.0E-02 

20 

12 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1.0E-07 

1 

The  choice  of  the  high  and  low  levels  was  based  upon  regional  values  for  massive 
Paleozoic  limestone  rocks  as  presented  by  Beck  and  others  (1988:  336-337).  The 
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approach  taken  was  to  use  as  wide  a  range  of  values  as  was  deemed  reasonable  for  the 
type  of  rocks  comprising  this  aquifer.  No  effort  was  made  to  limit  the  values  of  the 
hydraulic  conductivities  of  each  zone  to  any  particular  range,  although  such  delineation 
may  be  made,  if  desired.  The  single  porosity  value  used  represented  the  effective  porosity 
for  this  entire  fractured  carbonate  system. 

The  12  individual  experiments  were  run  through  SUTRA,  and  the  sum  of  the 
squared  differences  (SSE)  between  the  hydraulic  head  pressures  for  the  524  nodes  of  the 
run  and  the  calibration  target  set  was  computed.  To  determine  which  parameters  had  a 
significant  effect  on  the  output  response  (SSE)  of  the  model,  a  normal  probability  plot 
(Figure  4.4)  was  made  (Box  and  Draper,  1987: 128-134).  Using  this  method,  those 
factors  which  do  not  significantly  contribute  will  fall  toward  the  center  of  the  plot  and  lie 
along  a  line,  thus  indicating  they  do  not  appear  to  stand  above  the  normal  noise  of  the 
data.  In  figure  4.4,  Units  2,  9,  10,  and  possibly  6  appear  to  fall  off  of  the  line  near  its 
ends;  indicating  only  these  four  parameters  had  a  significant  effect  on  the  model. 

In  addition  to  this  graphical  technique,  a  step-wise  linear  regression  was  performed 
on  the  data.  Using  an  alpha  of  0.075  to  enter  and  exit  the  linear  model,  the  step-wise 
regression  procedure  also  selected  Units  2,  6,  9,  and  10  to  enter  into  the  model. 

This  screening  experiment  determined  that  only  four  parameters  needed  to  be 
varied  to  calibrate  the  finite-element  model.  Efforts  to  adjust  any  of  the  other  parameters 
would  have  minimal  effect,  as  their  contribution  is  no  greater  than  the  noise  contained 
within  the  data.  Therefore,  in  applying  the  RSM  to  the  calibration,  only  the  hydraulic 
conductivities  for  Units  2,  6,  9,  and  10  were  adjusted.  The  settings  for  all  the  other 
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parameters  were  fixed  to  values  at  the  center  of  the  Plackett-Burman  design.  Porosity 
was  fixed  at  10%  and  the  hydraulic  conductivities  for  Units  1,  3, 4,  5,  7,  and  8  were  fixed 
at  5.0  x  10'3  m/min.  It  should  be  emphasized,  however,  that  this  analysis  has  shown  that 
the  settings  for  these  units  are  not  significant,  so  they  may  be  set  to  any  value  in  the 
parameter  range. 

In  closing  this  discussion,  a  risk  associated  with  screening  designs  should  be 
highlighted.  Screening  designs  generally  provide  the  ability  to  estimate  only  linear  effects 
and  may  not  identity  certain  types  of  higher  order  effects.  For  example,  a  factor  that 
imparts  a  significant  quadratic  effect  on  the  response  could  potentially  be  included  in  the 
set  of  factors  having  no  appreciable  effect  on  it.  If  the  observed  responses  at  the  high  and 
low  levels  of  this  factor  were  approximately  equal,  the  estimated  linear  effect  would  be 
negligible.  The  important  quadratic  effect  would  have  been  overlooked,  as  a  design  with 
only  two  levels  cannot  estimate  a  quadratic  effect.  Thus,  it  is  possible  that  a  screening 
experiment  may  unknowingly  delete  one  or  more  variables  with  significant  higher  order 
effects  from  further  consideration. 

If  response  surface  methods  are  not  able  to  calibrate  the  finite-element  model 
sufficiently  well,  the  investigator  may  need  to  reexamine  or  reaccomplish  the  screening 
experiment  to  determine  if  significant  factors  were  inadvertently  omitted.  In  this  study,  it 
is  believed  that  no  influential  factors  were  incorrectly  screened  from  the  calibration  effort. 
Thus,  by  adjusting  only  the  four  hydraulic  conductivities  for  Units  2,  6,  9,  and  10;  a 
calibrated  model  should  be  achieved. 
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Figure  4.4  SSE  Normal  Probability  Plot 
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2666.63 

9 

0.7727 
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8 

0.6818 
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0.23 
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6 
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0.00 

Porosity 

35.18 
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5 

0.4091 

-0.23 

Unit  1 

28.08 
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4 

0.3182 

-0.4/ 

Unit  3 

-49.85 

-99.71 

3 

0.2273 

-0.74 

Unit  7 

-286.32 

-572.64 

2 

0.1364 

-1.09 

Unit  6 

-690.07 

-1380.14 

1 

0.0455 

-1.69 

SSE  Normal  Probability  Plot 


First-Order  Designs 

Once  the  significant  parameters  were  determined,  first-order  designed  experiments 
were  used  to  provide  local  approximations  of  the  response  surface.  Four  designs  (A,  B, 

C,  and  D)  were  used  to  fully  exploit  the  first-order  strategy. 

Based  upon  the  screening  experiments,  only  four  hydraulic  conductivity 
parameters  needed  adjusting  (the  other  six  hydraulic  conductivities  were  fixed  at  5.0  x  10"3 
m/min  and  effective  porosity  at  10%).  Thus  a  24  (two-level,  four  factor)  factorial  design 
was  used  throughout  the  first-order  phase.  This  design  required  16  runs  and  generated  the 
SSE  responses  that  were  regressed  to  estimate  the  gradient  at  each  design  region.  Tables 
of  the  specific  settings  used  for  each  of  the  four  24  factorial  designs  and  their  subsequent 
steepest  descent  trials  are  contained  in  Appendix  A. 

Design  A  The  starting  point  for  the  first  first-order  design  was  based  upon  the 
approximate  center  point  of  the  Plackett-Burman  design.  The  starting  point  could  have 
used  the  settings  from  the  best  run  found  during  the  Plackett-Burman  screening 
experiment  (Run  12,  SSE  =  0. 14596).  This  run  gave  a  rather  small  SSE  and  provided  a 
good  starting  point  for  calibrating  the  finite-element  model.  However,  the  intent  of  this 
study  was  to  demonstrate  the  general  applicability  of  the  RSM  technique.  The  discovery 
of  “good”  parameter  settings  was  not  exploited  in  order  to  demonstrate  the  methodology 
in  the  absence  of  such  fortuitous  results. 

Although  the  Plackett-Burman  design  used  the  same  settings  for  all  the  hydraulic 
conductivities,  the  deeper  units  were  assigned  smaller  values  in  Design  A  since  the 
hydraulic  conductivities  usually  decrease  with  depth.  The  high  and  low  settings  for  these 
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experimental  runs  were  based  on  ±  80%  of  the  center  point’s  value.  Table  4.3  presents 
the  four  hydraulic  conductivity  values  used  for  Design  A.  Table  4.4  presents  the  results  of 
regressing  the  16  SSE  responses  using  the  Regression  Analysis  Tool  from  Microsoft 
Excel.  The  Residual  statistic  was  partitioned  into  Quadratic  and  Other  Terms  as  presented 
in  Chapter  3,  “Single  Degree  of  Freedom  Test  For  Curvature.” 

This  first-order  model  had  high  explanatory  power  as  indicated  by  the  high  R 
Square  (0.9287)  and  F  value  (39.06).  However,  quadratic  terms  seemed  to  be  significant 
(as  seen  in  the  high  F  statistic  for  the  single  degree  test  for  curvature).  Because  of  the 
high  explanatory  power  of  this  first-order  model,  further  exploration  was  warranted  down 
the  path  of  steepest  descent  in  spite  of  the  evidence  of  curvature  suggested  by  the  single 
degree  of  freedom  test  for  curvature.  The  steepest  descent  gradient  was  determined  from 
the  negative  of  the  regression  coefficients  and  normalized  in  the  units  of  the  design  (Box 
and  Draper,  1987: 187).  Experiments  were  conducted  along  this  gradient  until  the  SSE 
response  failed  to  become  smaller.  Figure  4.5  illustrates  the  behavior  of  the  SSE  response 
along  this  gradient  determined  from  the  Design  A  experiments  (Note  the  curve  is 
interpolated  by  the  graphics  package;  only  the  points  represent  actual  data).  The  lowest 
SSE  response  occurred  at  14  unit  vector  lengths  from  the  center  of  Design  A.  The  SSE 
response  was  reduced  from  a  value  of  0.1573  to  0.1439.  The  point  14  unit  vector  lengths 
from  the  center  of  Design  A  served  as  the  center  of  the  next  design.  Design  B. 
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Table  4.3  Design  A  Experimental  Values 

Hydraulic  Low  Center  High 

Conductivity  Zone  (m/min)  (m/min)  (m/min) 


Unit  2 

1.0000  x  10'3 

5.0000  x  1  O'3 

9.0000  x  1  O’3 

Unit  6 

1.0000  x  lO'4 

5.0000  x  1  O’4 

9.0000  x  1  O'4 

Unit  9 

1.0000  x  10"6 

5.0000  x  1  O'6 

9.0000  x  lO45 

Unit  10 

1.0000  x  1  O'6 

5.0000  x  1  O'6 

9.0000  x  10‘6 

Table  4.4  Design  A  Regression  Analysis 


SUMMARY  OUTPUT 


Regression  Statistics 

Multiple  R 

0.9637 

R  Square 

0.9287 

Adjusted  R  Square 

0.9049 

Standard  Error 

1.347E-02 

Observations 

17 

ANOVA 


df 

SS 

MS 

F 

Significance 

F 

Regression 

4 

2.836E-02 

7.090E-03 

39.06459 

0.00000 

Residual 

12 

2.178E-03 

1.815E-04 

Quad.  Terms 

1 

1.532E-03 

1.532E-03 

26.08723 

0.00034 

Other  Terms 

11 

6.460E-04 

5.873E-05 

Total 

16 

3.054E-02 

Coefficients  Standard  t  Stat  P -value 

Error 


Intercept 

1.953E-01 

3.267E-03 

5.977E+01 

3.182E-16 

Unit  2 

2.919E-03 

3.368E-03 

8.668E-01 

4.031E-01 

Unit  6 

-4. 14  IE-02 

3.368E-03 

-1.229E+01 

3.691E-08 

Unit  9 

4.992E-03 

3.368E-03 

1.482E+00 

1.641E-01 

Unit  10 

-4.959E-03 

3.368E-03 

-1.472E+00 

1.666E-01 
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Designs  B,  C,  &  D  For  the  second,  third,  and  fourth,  first-order  designs,  the  high 
and  low  settings  were  fixed  at  ±  50%  of  the  center  point’s  value  in  an  attempt  to  reduce 
the  amount  of  curvature  in  the  resulting  empirical  models.  Tables  4.5  -  4.7  present  the 
hydraulic  conductivity  values  used  for  Designs  B,  C,  &  D.  Each  design  was  centered  on 
the  lowest  SSE  response  observed  during  the  steepest  descent  search  away  from  the 
previous  design. 


Table  4.5  Design  B  Experimental  Values 


Hydraulic 
Conductivity  Zone 

Low 

(m/min) 

Center 

(m/min) 

High 

(m/min) 

Unit  2 

5.5845  x  lO-4 

1.1 169  x  10‘3 

1.6754  x  10'3 

Unit  6 

3.0038  x  10‘3 

6.0076  x  10'3 

9.01 14  x  10’3 

Unit  9 

5.0000  x  10'8 

1.0000  x  10'7 

1.5000  x  1  O'7 

Unit  10 

5.7985  x  10"6 

1.1597  x  10‘5 

1.7396  x  10‘5 

Table  4.6  Design  C  Experimental  Values 

Hydraulic 

Low 

Center 

High 

Conductivity  Zone 

(m/min) 

(m/min) 

(m/min) 

Unit  2 

2.2201  x  10'2 

4.4402  x  1  O’2 

6.6603  x  10‘2 

Unit  6 

3.1991  x  10‘2 

6.3982  x  1  O'2 

9.5973  x  10‘2 

Unit  9 

1.2884  x  10’7 

2.5767  x  1  O'7 

3.8650  x  10'7 

Unit  10 

5.0000  x  10'8 

1.0000  x  1  O'7 

1.5000  x  10  -7 

Table  4.7  Design  D  Experimental  Values 

Hydraulic 

Low 

Center 

High 

Conductivity  Zone 

(m/min) 

(m/min) 

(m/min) 

Unit  2 

9.2195  x  10'3 

1.8439  x  1  O'2 

2.7658  x  10’2 

Unit  6 

6.2300  x  10‘2 

1.2460  x  10’1 

1.8690  x  10'1 

Unit  9 

7.7060  x  1  O'8 

1.5412  x  10'7 

2.3118  x  10’7 

Unit  10 

7.0100  x  10'8 

1.4020  x  10"7 

2.1030  x  10"7 
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For  each  design,  the  SSE  responses  were  regressed  on  the  coded  parameters,  the 
steepest  descent  gradient  was  determined,  and  experiments  were  run  along  this  gradient 
until  the  SSE  response  failed  to  become  smaller.  Figures  4.6  -  4.8  illustrate  the  behavior 
of  the  SSE  response  out  along  the  gradients  determined  from  Designs  B,  C,  &  D.  Both 
Designs  B  and  C  demonstrated  significant  improvement  (lowering  of  the  SSE  response), 
and  the  single  degree  test  for  curvature  indicated  only  a  minor  amount  of  curvature.  Table 
4.8  summarizes  the  results  from  all  of  the  first-order  designs. 


Table  4.8  Summary  of  First-Order  Designs 


Design 

R  Square 

Quadratic 
F- value 

Quadratic 

p-value 

Distance  From 
Design  Center 

Lowest  SSE  Response 
Along  Gradient 

A 

0.9287 

26.087 

0.00034 

14 

0.14394 

B 

0.8070 

4.127 

0.06707 

80 

0.11730 

C 

0.9216 

4.152 

0.06637 

2.25 

0.10318 

D 

0.8675 

94.156 

0.00000 

0.50 

0.10099 

Design  D  did  not  show  significant  improvement  over  Design  C,  the  R  Square  of 
the  fitted  first-order  model  fell,  a  very  significant  amount  of  curvature  existed,  and  the 
local  minimum  is  found  well  within  the  design  region.  These  observations  indicated  the 
first-order  approach  had  reached  its  limits  of  applicability  and  a  second-order  design 
appeared  justified  to  further  evaluate  the  response  surface. 

The  first-order  strategy  yielded  very  good  results;  improving  the  SSE  response 
from  0. 15732  (at  the  center  of  Design  A)  to  0. 10099  (after  following  the  gradient  from 
Design  D).  The  calibrated  model  0.50  units  from  the  center  of  Design  D  had  a  mean 
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absolute  error  (MAE)  of  only  7. 1505  mm  at  each  of  the  524  nodes.  In  fact,  the  calibrated 
model  obtained  after  searching  down  the  gradient  from  Design  C  would  probably  have 
been  good  enough,  in  most  cases,  to  use  in  groundwater  simulation  studies  without  further 
fine-tuning  of  the  input  parameters.  However,  to  fully  demonstrate  the  RSM  technique,  a 
second-order  model  was  evaluated  to  calibrate  the  SUTRA  model  even  closer  to  the 
calibration  target  data  set. 


✓ 
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B  Center  (Unit  Vector  Lengths) 


Figure  4.7  SSE  Response 
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Figure  4.8  SSE  Response  Al 
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Second-Order  Designs 

To  evaluate  the  SSE  response  surface  further,  a  second-order  model  was  fit  to  the 
data.  The  24  factorial  design,  used  in  Design  D,  was  augmented  with  8  axial  points  and  a 


center  point,  to  form  a  central  composite  design  (CCD)  experiment.  The  axial  length  (a) 
for  these  points  should  be  equal  to  2  in  order  to  have  a  rotatable  CCD.  However,  when 
a  =  2,  the  actual  parameter  values  became  zero.  To  avoid  this  problem,  an  a  of  1 .99  was 
used  in  the  first  CCD,  named  CCD  A.  Other  settings  of  a  were  used  to  investigate  the 
effect  this  setting  had  on  the  results  of  the  CCD  experiment.  Although  in  this  study  three 


different  values  of  a  were  used,  only  CCD  C  was  really  necessary.  Once  the  additional 
axial  points  were  run  through  SUTRA,  the  SSE  responses  were  regressed  to  a  second- 
order  model  using  the  SAS  statistical  analysis  computer  package.  SAS  also  presented  a 
ridge  analysis  of  the  fitted  surface,  suggesting  a  path  of  lowest  responses.  Table  4.9 
summarizes  the  results  of  the  three  CCDs. 

Table  4.9  Summary  of  Second-Order  Designs 


CCD  Name 

a  Value 

R  Square 

Lowest  Min.  Ridge 
SSE  Response 

CCD  A 

1.99 

0.4730 

0.10171 

CCD  B 

1.50 

0.9806 

0.09918 

CCD  C 

1.00 

0.9988 

0.09820 

CCD  C,  with  a  -  1.00,  had  the  best  fit  of  the  second-order  surfaces  with  an  R 
Square  of  0.9988.  By  following  the  suggested  minimum  ridge  away  from  the  center  of 
CCD  C  (based  upon  the  fitted  second-order  surface),  the  lowest  SSE  response  in  this 
study  was  found,  0.0982.  Figures  4.9  -4.11  show  the  SSE  responses  as  the  ridge 
minimum  is  followed  away  from  each  CCD. 
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Summary  and  Calibration  Precision 

Figure  4. 12  presents  the  lowest  SSE  responses  observed  from  (1)  the  Plackett- 
Burman  screening  experiment;  (2)  each  gradient  search  for  each  of  the  four  first-order 
models;  (3)  the  ridge  analysis  of  the  second-order  model,  CCD  C.  This  figure 
demonstrates  the  improvement  in  the  calibration  of  the  SUTRA  model  with  each  new 
designed  set  of  experiments.  Table  4. 10  summarizes  the  number  of  experimental  runs 
expended  to  achieve  the  various  levels  of  calibration.  As  discussed  earlier.  Design  A 
intentionally  was  not  set  at  the  best  response  observed  in  the  screening  design.  Selecting 
the  center  of  the  screening  design,  with  its  correspondingly  higher  SSE,  provided  a  greater 
opportunity  to  evaluate  the  appropriateness  of  both  the  SSE  response  and  the  RSM 
technique  for  calibrating  a  groundwater  flow  model.  Depending  on  the  level  of  precision 
desired,  the  calibration  effort  could  have  required  as  little  as  12  runs  or  as  many  as  146. 


Table  4.10  Number  of  Experimental  Runs 

Experimental  Design  Number  of  Runs 


Screening  Design: 

Plackett-Burman  Design  12 

First-Order  Designs: 

Design  A  16 

Gradient  Search  12 

Design  B  16 

Gradient  Search  30 

Design  C  16 

Gradient  Search  9 

Design  D  16 

Gradient  Search  5 

Second-Order  Design: 

CCD  C  8 

Min.  Ridge  Run  6 

TOTAL  RUNS  146 
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Figure  4.13  plots  the  best  computed  hydraulic  heads  (based  on  the  minimum  ridge 
run  1.75  units  away  from  CCD  C)  versus  the  calibration  target  set  of  heads  that  were 
being  matched.  The  data  plot  along  a  45°  line,  indicating  the  excellent  calibration  of  the 
SUTRA  model.  Table  4. 1 1  lists  other  measures  of  the  error  between  the  best  SUTRA  run 
and  the  calibration  target  set,  indicating  the  close  match  achieved. 


Table  4.11  Measures  of  Calibration  Precision 


Measure  of  Error 

Value  ( units) 

Sum  of  Squared  Error  (SSE) 

0.09819  m2 

Root  Mean  Squared  Error  (RMS) 

0.01369  m 

Mean  Absolute  Error  (MAE) 

0.00698  m 

Mean  Error  (ME) 

0.00148  m 

Maximum  AE 

0.07573  m 

Median  AE 

0.00197  m 

Minimum  AE 

0.00000  m 

Figures  4. 14  and  4. 15  plot  the  error  between  the  heads  computed  at  the  524  nodes  and  the 
horizontal  and  vertical  positions  of  the  nodes.  These  plots  show  the  errors  are 
symmetrically  distributed  over  almost  all  of  the  nodes,  indicating  the  SUTRA  model  has 
been  evenly  calibrated,  without  concentrating  its  errors  in  any  particular  area  of  the 
aquifer.  By  empirically  testing  the  SUTRA  model  and  systematically  evaluating  the 
response  surface  using  the  response  surface  methodology,  a  finely  calibrated  groundwater 
model  was  achieved. 

Although  the  RSM  technique  produced  a  very  well  calibrated  model,  the  final 
values  of  the  hydraulic  conductivities  found  do  not  match  the  values  used  in  the  Smith- 
Ritzi  model  very  closely.  Table  4.12  compares  the  settings  of  the  parameters  used  in  the 
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Smith-Ritzi  model  with  the  values  used  to  produce  the  calibrated  model  of  this  thesis. 
Recall  that  as  a  result  of  the  screening  experiment,  only  the  hydraulic  conductivities  of 
Units  2,  6,  9,  and  10  were  adjusted  using  RSM.  The  remaining  parameters  were  fixed  at 
notional  values.  The  differences  between  the  target  and  calibrated  parameter  settings 
illustrate  the  nonuniqueness  property  typical  of  many  inverse  problems.  However,  the 
goal  of  this  investigation  was  to  calibrate  a  groundwater  model  so  that  its  computed 
hydraulic  heads  closely  matched  the  calibration  target  data  set,  not  necessarily  the  true 
parameters  of  the  groundwater  model.  Clearly,  the  RSM  technique  has  located  a  local 
optimum  on  the  SSE  surface  that  results  in  a  finely  calibrated  groundwater  model  without 
having  matched  the  parameters  used  in  the  calibration  target  data  set. 


Table  4.12  Parameter  Settings:  Smith-Ritzi  Model  vs  Calibrated  Model 


SUTRA  Groundwater 

Model  Parameter 

Value  in 

Smith-Ritzi  Model 

Value  in 

Calibrated  Model 

Porosity 

11% 

10% 

Hydraulic  Conductivity  Unit  1 

2.000  x  10°  m/min 

5.000  x  10'3  m/min 

Hydraulic  Conductivity  Unit  2 

1.000  x  10°  m/min 

9.843  x  1 0'3  m/min 

Hydraulic  Conductivity  Unit  3 

1.000  x  10*4  m/min 

5.000  x  10'3  m/min 

Hydraulic  Conductivity  Unit  4 

3.000  x  10'5  m/min 

5.000  x  10'3  m/min 

Hydraulic  Conductivity  Unit  5 

2.000  x  10"3  m/min 

5.000  x  10'3  m/min 

Hydraulic  Conductivity  Unit  6 

3.000  x  10"4  m/min 

1.009  x  1 0'1  m/min 

Hydraulic  Conductivity  Unit  7 

6.000  x  10-4  m/min 

5.000  x  10'J  m/min 

Hydraulic  Conductivity  Unit  8 

6.000  x  10-6  m/min 

5.000  x  10"3  m/min 

Hydraulic  Conductivity  Unit  9 

4.000  x  10"6  m/min 

4.385  x  l O'8  m/min 

Hydraulic  Conductivity  Unit  10 

1.000  x  10'5  m/min 

1.400  x  10'7  m/min 
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SSE  Response 
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Figure  4.12  Best  SSE  Response 
From  Each  Designed 
Experiment 


Originating  Design 
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Figure  4.13  Computed  vs  Calibration  Target  Heads 
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Figure  4.15 
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V.  Conclusions  and  Recommendations 


Conclusions 

Response  surface  methodology  (RSM)  can  be  used  to  estimate  groundwater  model 
flow  parameters.  RSM  provides  an  efficient  and  effective  process  to  determine  those 
parameters  which  have  a  significant  effect  on  the  groundwater  model.  RSM  also  provides 
the  systematic  methods  to  adjust  these  parameters  to  quickly  arrive  at  a  calibrated  model. 

An  existing  two-dimensional,  steady-state,  saturated  flow  model  was  calibrated 
using  RSM.  A  screening  experiment,  based  upon  a  Plackett-Burman  design,  provided  an 
efficient  method  to  determine  the  most  important  groundwater  model  parameters  to 
adjust.  This  screening  experiment  gave  information  that  is  usually  obtained  through 
sensitivity  analysis  in  traditional  calibration  studies.  By  the  use  of  the  normal  probability 
plot,  a  direct  indication  of  the  statistically  significant  parameters  was  made  apparent. 

Once  the  significant  parameters  were  determined,  the  calibration  of  the  groundwater 
model  was  made  by  adjusting  these  chosen  parameters. 

Response  surface  methods  provided  an  effective  means  to  adjust  the  parameter 
values  to  calibrate  the  groundwater  model  to  the  calibration  target  data  set.  First-order 
models  were  exploited  to  rapidly  adjust  the  parameters  to  produce  a  well  calibrated  finite- 
element  model.  A  second-order  model  was  applied  to  further  calibrate  the  finite-element 
model  even  closer  to  the  calibration  target  data  set. 
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The  SSE  proved  to  be  a  suitable  response  for  the  RSM  calibration  technique. 
However,  the  fairly  significant  single  degree  of  freedom  tests  for  curvature  indicated  that 
the  SSE  response  surface  is  at  least  quadratic,  even  at  distances  somewhat  removed  from 
a  local  optimum.  Due  to  the  quadratic  properties  exhibited  by  the  SSE  surface,  the  design 
regions  for  the  first-order  designs  should  be  sufficiently  “small”  to  ensure  adequate  first- 
order  models.  As  an  alternative  measure,  the  root  mean  square  (RMS)  or  square  root  of 
SSE  would  probably  produce  a  “flatter”  response  surface.  With  a  flatter  response  surface, 
the  first-order  designs  may  converge  more  efficiently  and  may  preclude  the  need  to 
develop  any  second-order  models,  at  all. 

This  investigation  was  conducted  with  a  minimum  of  assumptions  regarding  the 
initial  values  of  the  hydrogeologic  parameters  to  use  in  the  groundwater  model.  In  this 
manner,  RSM  was  being  tested  to  see  if  its  methods  would  calibrate  a  groundwater  model 
in  a  realistic  setting.  The  groundwater  model  calibrated  with  RSM  may  be  used  to 
estimate  contaminate  transport  or  evaluate  flow  patterns,  now  that  the  hydraulic  heads  are 
calibrated. 

Response  surface  methods  provide  a  systematic  technique  to  model  calibration 
versus  manual  trial-and-error  methods.  RSM  enables  the  hydrologist  to  adjust  all  the 
parameters  simultaneously  to  arrive  at  the  calibration  settings.  By  adjusting  the 
parameters  simultaneously,  possible  interaction  effects  of  one  parameter  upon  another  are 
taken  into  account  as  they  are  mutually  adjusted.  Manual  trial-and-error  methods  that 
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only  adjust  one  parameter  at  a  time  are  likely  to  miss  interaction  effects  (Montgomery, 
1991:487). 

RSM  also  provides  a  good  alternative  to  the  automated  techniques  that  attempt  to 
solve  the  inverse  problem  by  deriving  the  shape  of  the  response  surface.  RSM  does  not 
need  any  special  computer  code  to  explicitly  derive  the  response  surface.  Using  any 
readily  available  spreadsheet  program,  the  local  nature  of  the  response  surface  may  be 
elucidated  through  first-order  linear  models.  Second-order  models  are  handled  most 
conveniently  with  the  S  AS  statistical  package,  but  may  not  be  required  since,  in  some 
cases,  the  groundwater  model  may  be  adequately  calibrated  using  only  first-order  designs. 

In  summary,  RSM  has  been  demonstrated  to  be  a  viable  method  of  calibrating 
groundwater  flow  models.  Its  methods  can  be  readily  applied  in  the  calibration  process  to 
both  choose  the  relevant  parameters  to  adjust  and  to  determine  at  what  level  they  should 
be  set  to  achieve  a  calibrated  model. 

Recommendations 

The  following  are  some  recommendations  for  extending  this  research: 

1 .  Investigate  how  well  the  calibration  proceeds  in  a  data-poor  situation.  Rather  than 
using  a  calibration  target  data  set  of  524  nodes,  a  calibration  target  data  set  with  about 
12-18  observation  points  (wells),  or  less,  would  be  more  representative  of  a  “typical” 
water  well  field  (Carrera  and  Neuman,  1986c).  Kriging  techniques  may  also  be 
investigated  to  enrich  the  limited  data. 
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2.  The  initial  starting  values  of  the  groundwater  model’s  parameters  will  determine  how 
closely  the  groundwater  model  may  be  calibrated  to  the  calibration  target  data  set. 

The  inverse  problem  is  solved  using  RSM  to  locate  a  local  minimum.  The  influence  of 
different  starting  values  on  the  minimum  actually  found  should  be  evaluated. 

3.  Instead  of  using  the  SSE  response  surface,  a  possibly  “flatter”  surface  may  be  realized 
by  using  RMS  or  the  square  root  of  SSE.  If  this  surface  is  flatter,  first-order  strategies 
may  be  exploited  with  even  greater  efficiency. 

4.  If  the  groundwater  flow  model  is  to  be  used  in  contaminant  transport  studies,  it  should 
be  calibrated  for  solute  transport,  too.  Joint  optimization  can  be  handled  by  RSM 
techniques  and  should  be  investigated  for  a  groundwater  modelling  problem  with 
simultaneous  flow/transport  calibration  requirements. 
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Appendix  A 


Experimental  Design  Settings  and  Responses 

This  appendix  contains  the  Microsoft  Excel  worksheets  used  to  record  the 
parameter  settings  for  the  experimental  runs.  The  experimental  settings  and  the  responses 
for  each  run  are  listed  for  the  Plackett-Burman,  four  first-order,  and  three  second-order 
designs.  Also  included  are  the  gradient  or  ridge  analysis  searches  away  from  the  designed 
regions. 
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Appendix  B 


SUTRA  FORTRAN  Post-Processor  File 

This  appendix  contains  the  FORTRAN  source  code  used  to  automatically  compute 
the  error  statistics  that  measured  the  differences  between  the  SUTRA  output  and  the 
Smith-Ritzi  calibration  target  data  set.  This  code,  named  POST.FOR,  was  automatically 
called  by  the  VAX  command  file  in  Appendix  C  to  produce  the  SUTRA  output  report 
after  each  model  execution. 
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PROGRAM  POST 
C 

C  THIS  PROGRAM  PERFORMS  POST  PROCESSING  OF  THE  SUTRA  OUTPUT 

C  TO  DETERMINE  VARIOUS  ERROR  STATISTICS. 

C 


C 


DIMENSION  HBASE (600) , HNEW (600) 

REAL  SSE, RMS, MAE, ME 

OPEN ( 8 ,  FILE= ' HBASE . dat ' , STATUS= ’ UNKNOWN ’ , FORM= ' FORMATTED ’ ) 
OPEN (9,  FI LE=' node 8 3. dat' , STATUS= ' UNKNOWN ' , FORM= ’ FORMATTED ' ) 
OPEN (10,  FILE= ’ input . dat ' , STATUS= ' UNKNOWN ’ , FORM= ’ FORMATTED ' ) 
OPEN (11,  FILE= ' final . rpt ’ , STATUS= ' UNKNOWN ' , FORM= ’ FORMATTED ' ) 


WRITE (11, 90)  '  ' 
WRITE (11, 90)  ' 

ssss 

UU 

uu 

TTTTTT 

RRRRR 

AA' 

WRITE (11,90)  ' 

ss  s 

uu 

uu 

T  TT  T 

RR  RR 

AAAA' 

WRITE (11, 90)  ’ 

ssss 

uu 

uu 

TT 

RRRRR 

AA  AA' 

WRITE (11, 90)  ' 

ss 

uu 

uu 

TT 

RR  R 

AAAAAA' 

WRITE ( 11, 90)  ’ 

ss  ss 

uu 

uu 

TT 

RR  RR 

AA  AA’ 

WRITE (11,90)  ' 

ssss 

uuuu 

TT 

RR  RR 

AA  AA’ 

WRITE (11,90)  ’  ' 

+  *  ************  *  *  < 
WRITE (11,  90) 

WRITE (11, 90) 

WRITE (11, 90) 

WRITE (11,  90) 


lmeasurements  are 


WRITE (11,  90) 

1.  Smith' 's  thesis 
WRITE (11, 90) 

WRITE (11, 90) 

1 - 

WRITE (11, 90) 

WRITE (11, 90) 

1 - 

WRITE (11, 90) 

WRITE (11, 90) 

90  FORMAT  (IX, A) 

WRITE ( 11, 90)  ' 


SUBSURFACE  FLOW  SIMULATION  MODEL' 


Output  report  for  R.  M.  Cotman’’s  MS  thesis.  Error 


compared  to  the  steady  state  heads  presented  in  R.  J 


INPUT ' 


Input  SUTRA  data  file:  filename.D5' 


READ (10,*)  ITEMP 
WRITE (11, 110)  ITEMP 

110  FORMAT (IX, 'Porosity  (Percent):  ’,13) 

WRITE (11,90)  '  ' 

DO  30  1=1,10 

READ (10,*) TEMP 
WRITE ( 11, 100 ) I, TEMP 

100  FORMAT (IX, 'Hydraulic  conductivity  for  Unit', 13,'  (m/min) :  ',E10 

1.4) 

WRITE (11,90)  '  ' 

30  CONTINUE 

DO  10  1=1,524 

READ ( 8 , * ) X, Y, HBASE ( I ) 

READ ( 9, * ) X, Y, HNEW ( I ) 

10  CONTINUE 
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c 

WRITE (11, 90)  '  ' 

WRITE (11, 90)  ' - 

1 - . 

WRITE (11, 90)  '  OUTPUT' 

WRITE (11,90)  ' - 

1 - ' 


WRITE (11,90)  1  ' 

C  -  Sum  of  Squared  Error  (SSE)  Computation 

SSE=0 . 0 
DO  20  1=1,524 

SSE=SSE+ (HBASE (I ) -HNEW (I ) )**2 
20  CONTINUE 

WRITE (11, 120) SSE 

120  FORMAT (IX, ’Sum  of  Squared  Error  (SSE):  ' ,E15.5) 

WRITE (11, 90)  '  ' 

C  -  Root  Mean  Squared  Error  (RMS)  Computation 

RMS=0 . 0 

RMS=SQRT (SSE/524) 


WRITE (11, 122) RMS 

122  FORMAT ( IX, ' Root  Mean  Squared  Error  (RMS ) : ' , E15 . 5 ) 
WRITE (11, 90 )  '  ' 

C  -  Mean  Absolute  Error  (MAE)  Computation 

MAE=0 . 0 
DO  40  1=1,524 

MAE=MAE+ABS ( HBASE ( I ) -HNEW (I ) ) 

40  CONTINUE 

MAE=MAE/524 
WRITE (11, 12 4) MAE 

124  FORMAT (IX, 'Mean  Absolute  Error  (MAE):  ’,E15.5) 

WRITE (11, 90 )  '  ' 

C  -  Mean  Error  (MAE)  Computation  — 

ME=0 . 0 

DO  50  1=1,524 

ME=ME+ (HBASE ( I ) -HNEW { I ) ) 

50  CONTINUE 

ME=ME/524 
WRITE (11, 12 6) ME 

126  FORMAT (IX, 'Mean  Error  (ME):  ’,E15.5) 

WRITE (11, 90 )  '  ' 

STOP 

END 
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Appendix  C 


VAX  Command  File 


This  appendix  contains  the  VAX  command  file  used  to  simplify  the  execution  of 
the  SUTRA  model.  The  command  file,  named  SUTRA.COM,  served  as  an  interactive 
interface  to  the  SUTRA  program.  To  invoke  the  command  file,  the  user  simply  entered 
“@SUTRA”  at  the  VAX  command  prompt.  Upon  execution,  the  command  file  would 
prompt  the  user  for  the  porosity  value  and  the  settings  of  the  ten  hydraulic  conductivities. 
Once  these  parameters  were  entered,  the  SUTRA  input  parameter  file  was  automatically 
created,  the  SUTRA  program  was  executed,  and  the  error  statistics  were  computed  using 
the  POST.FOR  program  contained  in  Appendix  B,  which  also  created  an  output  report. 
This  output  report  was  saved  as  a  file  in  the  current  directory  and  displayed  to  the  screen 
for  immediate  review.  A  sample  output  report  is  presented  as  the  last  page  of  this 
appendix.  It  presents  the  results  from  the  run  using  the  parameters  which  produced  the 
minimum  SSE  response. 
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*  " 
*  >' 
*  " 
*  •' 


$write  sys$output  "***********+****************************" 
$write  sys$output  "*  SUTRA  Interactive  Data  Input  Program  *" 

$write  sys$output  "*  For  Capt  R.  Cotman's  MS  Thesis 

$write  sys$output  "*  Answer  every  question  in  the  units 
$write  sys$output  "*  shown  in  the  (  ) ,  using  the  format 

$write  sys$output  "*  shown  by  the  [  ] . 

$write  sys$output  "****************************************" 
$inquire  pi  "Filename  of  the  input  file  you're  creating  [No 
extension] : " 

$input_file  =  pi 

$inquire  pi  "Value  for  Porosity  (percent)  [xx] " 

$porosity  =  pi 
$cnt  =  0 


open  /write  file  unit.tmp 
open  /write  filel  input.dat 
write  file  "$edit  template. d5" 
write  file  "AZ" 

write  file  "sub/po/ ' 'porosity'/w" 
write  filel  '"'porosity'" 


1 


"Hydraulic  Conductivity  for  unif'cnt'  (m/min)  [x.xxxxE- 
""  then  goto  finish 


$ 

$ 

$ 

$ 

$ 

$ 

$loop : 

$cnt  =  cnt  + 

$inquire  pi 
xx]  " 

$if  pi  . eqs. 

$string  =  pi 

$  write  filel  '"'string'" 

$  write  file  "sub/  unit '' cnt '/’’ string ' /w" 

$if  cnt  .eqs.  10  then  write  file  "sub/  unit '' cnt '/’’ string ' /w" 
$if  cnt  .eqs.  10  then  write  filel  ’"'string'" 

$if  cnt  .eqs.  10  then  goto  finish 

$goto  loop 

$finish: 

$  write  file  "exit" 

$close  file 
$close  filel 
$@unit.tmp 

$  rename  template.  d5  '  input__f  ile .  D5 
$  open  /write  file  unit.tmp 

$  write  file  "$edit  sutemp.fil" 


write  file  "AZ" 

write  file  "sub/INPUT/ ' ' input_file ' /w" 
’exit" 


$ 

$ 

$  write  file 

$close  file 
$@unit. tmp 

$rename  sutemp.fil  SUTRA. FIL 

$cls 

$cls 

$write  sys$output  "Please  wait  about  15  seconds,  while" 
$write  sys$output  "the  SUTRA  model  runs." 

$run  main 
$run  post 

$  open  /write  file  final. tmp 

$  write  file  "$edit  final. rpt" 

$  write  file  "AZ" 

$  write  file  "sub/filename/ ’’ input_file ’ /w" 

$  write  file  "exit" 
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$close  file 
$@final . tmp 

$rename  f inal . rpt  ' input_file . RPT 
$del/noconfirm  input.dat;* 

$del/noconfirxn  unit. tmp;* 

$del/noconfirm  final. rpt;* 

$del/noconfirm  final. tmp;* 

$pu  sutra. fil 
$ty  ' input_file . RPT 

$write  sys$output  "You're  SUTRA  model  has  run.  The  " 
$write  sys$output  "output  report  is  in  the  file  named:" 
$write  sys$output  " ' ' input_file ’ . RPT” 

$exit 
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**************************************************************************** 

SUBSURFACE  FLOW  SIMULATION  MODEL 

Output  report  for  R.  M.  Cotman's  MS  thesis.  Error  measurements  are 
compared  to  the  steady  state  heads  presented  in  R.  J.  Smith's  thesis. 


INPUT 


Input  SUTRA  data  file:  VB175.D5 
Porosity  (Percent) :  10 


Hydraulic 

conductivity 

for 

Unit 

1 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

2 

(m/min) : 

0 . 9843E-02 

Hydraulic 

conductivity 

for 

Unit 

3 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

4 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

5 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

6 

(m/min)  : 

0. 1009E+00 

Hydraulic 

conductivity 

for 

Unit 

7 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

8 

(m/min) : 

0 . 5000E-02 

Hydraulic 

conductivity 

for 

Unit 

9 

(m/min) : 

0.4385E-07 

Hydraulic 

conductivity 

for 

Unit 

10 

(m/min)  : 

0. 1400E-06 

OUTPUT 


Sum  of  Squared  Error  (SSE) : 
Root  Mean  Squared  Error  (RMS) : 
Mean  Absolute  Error  (MAE) : 


0 . 98201E-01 
0 . 13690E-01 
0 . 69783E-02 


Mean  Error  (ME) : 


0. 14810E-02 
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